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ABSTRACT 

Using both numerical and analytical approaches, we demonstrate the existence of an 
effective power-law relation L oc vnP between the mean Lyapunov exponent L of stellar 
orbits chaotically scattered by a supermassive black hole in the center of a galaxy and 
the mass parameter m, i.e. ratio of the mass of the black hole over the mass of the 
galaxy. The exponent p is found numerically to obtain values in the range p ~ 0.3- 
0.5. We propose a theoretical interpretation of these exponents, based on estimates of 
local ‘stretching numbers’, i.e. local Lyapunov exponents at successive transits of the 
orbits through the black hole’s sphere of influence. We thus predict p = 2/3 — q with 
q « 0.1-0.2. Our basic model refers to elliptical galaxy models with a central core. 
However, we hnd numerically that an effective power law scaling of L with m holds 
also in models with central cusp, beyond a mass scale up to which chaos is dominated 
by the influence of the cusp itself. We finally show numerically that an analogous law 
exists also in disc galaxies with rotating bars. In the latter case, chaotic scattering by 
the black hole affects mainly populations of thick tube-like orbits surrounding some 
low-order branches of the xi family of periodic orbits, as well as its bifurcations at 
low-order resonances, mainly the Inner Lindbland resonance and the 4/1 resonance. 
Implications of the correlations between L and m to determining the rate of secular 
evolution of galaxies are discussed. 
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1 INTRODUCTION 


There is by now overwhelming evidence that black holes with masses ranging from 10® to lO^^Mfn (see review Ferrarese & FordI 

(2005 

11 exist in the center of most galaxies Kormendv et al.l ( 

I995II. Gebhardt et al.l (Il996ll. Faber et al. (Il997ll. Kormendv et al. 

(1997 

1. Kormendv et al. ( 

19981. Ivan der Marel et al.l (Il997 

1. Gebhardt et al.l (2000ll. Giiltekin et al.l (l2009al. Giiltekin et al. 

(2009b 1. Kormendv fc He 

(I2OI3I1. McConell & Ma (20131 

The presence of a supermassive black hole leads to a number 


of consequences for dynamics, whose study has b een a subject of ex t ended research in the last three d ecades (indicative 


refere n ces closely related to our present work are GerhaxW&^inn^j[l9^ /j_ Merritt&FYidman _([]h9€ )j_ Mgmtyfc Wnuril 


1 1996tl . Fridman &: Merritt ( 1997 1. Kandrup &: Sideri J ( 2002h . Kalapotharakos et al.l ( 2004 1. Kalapotharakos fc VoglisI ( 200^ 


see Merritt (1999) or (2013) for complete reviews of the subject.) 

Inter alia, the growth of a central black hole provides a mechanism driving secular evolution in galaxies. The creation 


cal galaxies JMerritt fc Ouinlan d 

1998t), 

Hollev-Bockelmann et al.l (l2001^. Hollev-Bockelmann et al.l ( 

2002h. Contopoulos et al. 

(I2OO2II. Kalapotharakos fc Voglis 

(I2OO5 

1. Muzzio et al.l (2005h. .lesseit et al.l (I2OO5II. Muzziol J200fitl. KalaootharakosI (I2OO8II. 

Valluri et al. (2O10ll. Vasiliev & Athanassoulal (2OI2III. The secular evolution causes two main effects: il the densitv profile 
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of the system by N-body simulations, it is found that the newly formed chaotic orbits start inducing a gradual change in the 
distribution of matter, i.e., the shape of the system, which becomes more spherical given that the distribution of chaotic orbits 
is more isotropic. The slow change in shape, in turn, causes an adiabatic change in the potential, thus affecting the phase 
space structure, in particular at energies for which the phase space was initially (before the insertion of the central mass) 
occupied mainly by box orbits passing arbitrarily close to the center. The main change of the phase space structure regards 


the volume increase of the domain corresponding to regular short axis tube (SAT) orbits fsee lBinnev fc Tremaina (120081 1 for 


definition). As the volume of the SAT domain in creases, some chaotic orbits are adiabatically captured inside the boundary 


of the SAT domain IIKalapotharakos et al.l (120041 11. This capture then induces an additional change in the form of the system, 


enhancing the conversion of chaotic orbits into SAT orbits, etc. This cyclic process maintains secular evolution up to a point 
where the population of chaotic orbits decreases substantially. The systems evolved by this mechanism are closer to oblate. 
Furthermore, in the final stages the percentage of chaotic orbits becomes smaller even than the one in the original systems, 
before the insertion of the central mass. 

It should be noted that the degree up to which the transformation of a system from triaxial to axisymmetric proceeds 
depends on how many box orb its are transformed to chaotically scattered orbits due to the central mass. For example, in 
Hollev-Bockelmann et al. ( 2 OO 2 II the initial percentage of outer box orbits is such that an adia batic introduct i on of the central 


mass does not destroy triaxiality in the outer parts of their galaxy simulation models. Also, Valluri et al. ( 201fll 'l examined 


how stochastic (or ‘non-reversible’) can the whole process of secular evolution be characterized, by considering the secular 
evolution of dark matter haloes in various models of central mass concentrations. While their findings re-confirm the picture 
of (non-reversible) stochastic evolution in the case of a point-like central mass concentration (e.g. a black hole), they find that 
there is also a different, i.e., ‘regular’ (or reversible) type of evolution taking place in systems in which the central mass is 
spatially distributed (e.g. a galactic disc or bulge embedded in a triaxial halo). 

Hereafter, we focus on the mechanism of secular evolution induced by the c haotic scattering of cen trophylic orbits. In 
order to quantify the rate of secular evolution induced by the above mechanism, Kalanotharakosl ( 2008h introduced a novel 
quantity, measurable in N-body simulations, called the effective chaotic momentum 


C = 


^total 


( 1 ) 


In Eq.(ITJ, Lu, is the mean value of the Lyapunov Characteristic Number (LCN) of a sub-ensemble of chaotic orbits in the 
system after the insertion of the central mass. This is defined by the orbits belonging to a percentage window ±0.3 around 
the characteristic value where the distribution P{\ogLj) of the logarithms of the Lyapunov exponents Lj of all the particles 
in chaotic orbits presents its global maximum. Considering a percentage window is necessary since, as we will see also below, 
the distribution P(logLj) is two-peaked, while the secular evolution is caused mainly by orbits forming the higher of the two 
peaks of the distribution. As a rule, these are centrophilic orbits passing arbitrarily closely to the central mass. On the other 
hand, AAu, is the total mass inside the same window, while Ntotai is the total mass of the N-body system considered. 

The use of the effective chaotic momentum £. allows one to quantify several phenomena related to the rate of secular 
evolution. A remarkable result is that there appears to be a global (i.e. the same in all simulations) theshold value Cth such 
that as a system undergoes secular ev olution, with a t i me-va rying value of C{t), the evolution becomes inefficient over a 
Hubble time when C{t) falls below Cth ( Kalanotharakosl ( 2008I H. 

From the two factors in the definition of C (Eq[T]), the percentage of chaotic orbits /Ntotai depends on the specific 

system s tudied, i.e. on the init i al dis tribu tion function that dete rmines the initial conditions of the simulation. However, as 
noted in Kalapotharakos et al. (2004) and Kalanotharakosl ( 2008ll . the distribution of Lyapunov exponents, and in particular 
the value of Lu, is found numerically to be not very sensitive to the choice of initial distribution function in the simulations. 
Thus, simulations representing systems with different profiles and triaxiality were found to exhibit different percentages of 
chaotic orbits, but similar levels of L™, the latter found, instead, to be correlated with the value of the mass ratio of the 
central mass over the mass of the hosting system. Hereafter, this ratio is simply referred to as the ‘central mass parameter’m. 
These findings indicate that, from the two factors entering the computation of the rate of secular evolution via the ‘effective 
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chaotic momentum’, one (percentage AN^/Ntotai of chaotic centrophylic orbits) depends on the self-consistent distribution 
function of the system considered, while the other (mean Lyapunov exponent L™) depends strongly only on the value of the 
mass parameter m. 


In the present paper we focus on this latter dependence, and seek to model its dynamical origin. We note that a dependence 
of t he mean Lyapunov e xpon e nt L of the centrophilic o rbits on m is a resu l t der ived also in studies u s ing fi xed potentials 
(e.g. Gerhard fc Binnev ( 1985 1. Merritt fc Valluri ( 1996I L Kandrup fc Siderli ( 2 OO 2 I I'). In Kalapotharakos ( 20081 ') . on the other 
hand, this dependence was explicitly determined by the orbital data of the particles in N-body simulations yielding the value 
of Lu, which is a good measure of L. A remarkably good power-law dependence was found, namely 


L oc mF, p Ks 0.5 


( 2 ) 


The proximity of p to 0.5 was also found in models of simple galactic potentials KalapotharakosI ( 2008l l. As shown 
2, somewhat smaller value s, aroun d p ~ 0.3, are deduced by a careful a posteriori analysis of the numerical results 


Merritt fc Valhiril ( lOOfih and in Kandrup fc Sideris ()200 


in section 
presented 


In the present paper we first reconfirm numerically the power-law m in further fixed-potential computations. Then, we 
develop a theoretical modelling allowing to interpret the origin of this power-law. We also justify why the exponent p obtains 
values in the observed range. In particular, our theory yields an exponent p = 2(3 — q, where q ~ 0.1-0.2. 


Our theoretical modeling stems from considering the dynamics of chaotic centrophilic orbits which undergo ‘transits’, 
i.e. they spend par t of th eir time inside and another part outside the so-called ‘sphere of influence’ of the central mass 
( Gerhard fc Binnev ( 1985l lf . One may note here that transi ting is a necessary c ondit ion for chaos, since orbits lying entirely 
within the sphere of influence (i.e. the so-called ‘pyramids’ ( Merritt fc Vasili^ (l201ll ')L obey three quasi-integrals of motion 
which are deformations of the Keplerian integrals derived using perturbation theory (here, the perturba tion consists of the 
triaxia l galactic potential, which perturbs the otherwise Keplerian dynamics very close to the black hole, see lMerritt fc Vasiliev 
( 201lh l. On reverse, as shown below, for transiting orbits one can determine the so-called stretching number (i.e. local Lyapunov 
number) yielding the local rate of deviation of two orbits with nearby initial conditions across one transit. The Lyapunov 
exponent, modeled as the sum of many such stretching numbers, turns then to be positive, i.e. the orbits are chaotic. We 
provide various types of evidence for the validity of this approximation, which allows to predict a (positive) mean value for the 
Lyapunov exponent as a function of the central mass parameter m. It is remarkable, in this respect, that during a transit the 
motion can be characterized as nearly Keplerian, while far from transits the motion is box-like and obeys three approximate 
integrals. Thus, both ‘piecewise’ motions can be analyzed by nearly-integrable approximations. Nevertheless, their union yields 
chaos. 


In our basic modeling we use a simple galactic model with a harmonic core. This ensures a priori the presence of many box 



chaos for centrophylic orbits. Thus, they significantly affect the value of L even without the presence of a central black hole. 
We will show, however, by numerical tests, that the presence of the central cusp introduces a critical mass parameter scale nic, 
depending on the value of 7, such that, for m > nic the chaotic scattering is dominated by the black hole, while for m < 
it is dominated by the central cusp. As shown in section 4, we then essentially recover an effective power law L oc nnF for 
m > rUc- Finally, we present numerical evidence that an effective power-law of the same form applies in rotating disc-barred 
galaxies with central black holes. In this case, the relevant mass parameter m corresponds to the ratio of the mass of the 
black hole over the total mass of the bar. In summary, although our theoretical interpretation for the origin of the effective 
law L oc is strictly valid in a quite simplified (and rather unrealistic) galactic model, our numerical evidence is that such 
a law appears generically in models of increasing complexity (and astrophysical interest). 

The structure of this paper is as follows: section 2 presents further numerical results about the empirical relation L oc 
using a simple galactic-type potential to which we add the potential of the central mass. These results are used in order to 
probe numerically some aspects of subsequent modeling. Our theoretical modeling itself is presented in Section 3. Section 4 
presents numerical results in models with central cusps and discusses the extent and limits of validity of previous results on 
the correlation between L and m in such models. Section 5 deals with the L oc relation in barred disc galaxy models, 
discussing both its applicability and origin, despite the non-existence, in such systems, of box-like orbits. Section 6 summarizes 
the main conclusions of the present study. 


























































































4 N. Delis, C. Efthymiopoulos and C. Kalapotharakos 


2 NUMERICAL RESULTS 

2.1 Hamiltonian model and numerical integrations 


At first we study the relation between L and m in a simple Hamiltonian model that captures the main features of motion 
near the center of an elliptical galaxy with non-singular central force held, to which we add a Keplerian term corresponding 
to the central mass. The Hamiltonian is: 


( 3 ) 


( 4 ) 


H{x,y,z,pa,,Py,Pz) = Y + y + Y + 
where V = V{x, y, z) is the gravitational potential: 

V(x, y, z) = + ^(.xz + xy + y'^f - mj+ d'^ . 

The variables {x,y,z) are cartesian position coordinates, {px,Py,Pz) their conjugate momenta, and r = . The 

softening parameter d was added in the Keplerian potential in order to avoid large numerical errors when the orbits pass very 
close to the center. We selected a set of incommensurable frequencies lo\ = 1, liJ 2 = V^, cua = so as to avoid dealing with 
resonant orbits satisfying some low-order commensnrability. The anharmonicity parameter e was given values between 0.01 
and 0.5. Various tests of the robustness of our results against e are reported below. We also note that the quartic potential 
term was choosen so as to represent a generic form without particular symmetries, while in galaxies with one or more planes 
of symmetry the potential presents a corresponding even symmetry. 

The gravitational potential 0 corresponds to a galaxy model with a harmonic core. The far more realistic case in which 
a central cusp is present is examined in detail in section 4 below. Here, the choice of the potential ensures a priori the existence 
of many regular box orbits, when m = 0, which are transformed to chaotically scattered orbits when m ^ Q.\n the model 0, 
the force grows linearly with distance from the center, while the force from the ce ntral mass falls as an inve rse square law. 
The two forces become similar in magnitnde at distances comparable to the radius I Gerhard fc Binnev (198^)): 


1/3 

= m ' Tc 


( 5 ) 


The sphere r = Vm is called sphere of influence of the central mass. The parameter is of order unity in units in which the 
frequencies uJx, <^y, are of order unity (as in Eq.ljS])). Then, for the total mass of the galaxy one also has M ~ r® = 0(1). 
The periods of orbits reaching apocentric positions Va » Tm are of order T ~ 27r. 

The equations of motion resulting from the Hamiltonian 0 are: 

X ^ Px 
y = Py 
z ~ pz 

Px = —ojlx — 2e{xz + xy + y'^){z + y) + mx/{D + d^)^^^ 

Py = -uly - 2 £{xz + xy + y^){x+ 2y) + my/{D + d^f^^ 

Pz = —UJ3Z — 2e(x2 + xy + y^)x + mzHp -|- . (6) 

In Lyapunov exponent computations we also make use of the variational equations: 


5x 

11 



Sy 

11 



Sz 

II 



Spx 


ddV . 
dxdy^^ 

d^V . 

_ Q2, 

dxdz 

Spy 

ddV _ 
dydx^"^ 

- -TT^y- 
dy 

d^V 

dydz 

5pz 


d^V - 
dydz^^ 

- 

dz^ 


( 7 ) 


In numerical computations we solve together the equations and ©. We use a 7-8th order Runge-Kutta method with fixed 
time step At = 10“^. For a choice of softening d = 10“®, this time step ensures that energy is conserved to within an error 
between 10“^^ and 10“^° for all orbits. 

The Lyapunov characteristic number (LCN) is defined through the relation: 


LCN= hm iln(C(t)/e(0)) 

t—^oo t 


( 8 ) 
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Figure 1. (a)Time evolution of finite Lypunov numbers for the ensemble of orbits in the case m = 10“^, s = 0.1. (b) The detailed form 
of the curve x(i) for one chaotic orbit (initial conditions x = 0.064, y = —0.0625, .2 = —0.123, Vx = 0.064, Vy = 0.028, Vz = —0.06. After 
the time t = 10^ the variations of x(f) ^-re less than 2%. 


where ^(t) = In practical computations we make use of the finite 

time Lyapunov characteristic number: 

x(t) = iln(C(t)X(0)) . (9) 

We made computations for ensembles of n = 200 orbits for the mass parameter values m = 10“®, 3 x 10“®, 10“^^, 
3 X 10”“^, 10“®, 3 X 10“®, 10“^, as well as the values of the anharmonicity parameter e = 0.01, 0.1 and 0.5. The initial 
conditions of each ensemble are choosen as follows. For each initial condition, we choose first randomly a value of the energy 
in the range 0 ^ E ^ 0.2, with uniform distribution. Then, we produce an orbit of zero initial angular momentum, by setting 
Px = Py = Pz = 0 and by solving for r the equation E = V{r,9,<j)), where {r,6,(f>) are spherical coordinates corresponding 
to the cartesian point {x,y,z), and cos6', (f) are choosen randomly with a uniform distribution in the intervals [—1,1] and 
[0, 27r) respectively. The selected range of energies represents motions with apocentric distances of order r « r^. These are all 
centrophilic orbits with a zero mean value of either component of the angular momentum. 


2.2 Results 


Figure [T] shows the time evolution of the finite-time Lyapunov exponent x(f) for the whole ensemble of orbits in one of the 
numerical experiments where m = 10“® and e = 0.1. The integration was up to the time t = 10®. Almost all orbits in this 
ensemble are chaotic, as their Hnite-time Lyapunov exponents are stabilized to non zero values. Only a small subset of orbits 
exhibit a value of x(t) that keeps falling even at the time t = 10®, following the law x(t) ~ 1/t- This subset defines regular 
orbits. The detailed time behaviour of logx(f) for a typical chaotic orbit is shown in FiglTJr. A general remark is that for 
the entire set of parameter values used in our experiments, the large majority of orbits in our ensembles turn to be chaotic. 
The minimum percentage of chaotic orbits (67%) is observed in the experiment with the minimum values of m and e, i.e. 
m = 10~® and e = 0.0 1. The classification of orbits as ordered or chaotic is based on a ‘Fast Lyapunov Indicator’ criterion 
i Froeschle et al. I IOOtI IL namely on whether the length ^(t) of the deviation vector is smaller, or larger, respectively, than a 
threshold value set equal to ^th{t) = lOOt. 

Figure [2] shows the distribution of the quantity logx(f) for three different time snapshots [t = 10®, t = 10“^, t = 10®) for 
the ensembles of orbits in the experiments with central mass values m = 10 ~^, m = 10 ~® and m = 10 “^. In all cases we can 
see that the distribution logx(^) exhibits a main peak corresponding to the chaotic orbits, which is displaced to higher values 
of logx(i) as the central mass m increases. On the other hand, for the mass values m = 10 “’’ and m = 10 “® there appears a 
secondary peak in the distribution of logx(i), that corresponds to the small subset of regular orbits. The secondary peak is 
displaced to the left, as the quantity log x(f) for regular orbits decreases in time as — log t. For the largest central mass values 
(m = 10“^), however, we observe no secondary peak, i.e., all the orbits turn to be chaotic. 

Figure [3] shows the main result. From the histograms of Fig[2] the mean value L = x(t) is extracted and plotted against 
m at the snapshots t = 10®, t — 10"* and t = 10® for all the numerical experiments. The straight lines in the same plots 
(in logarithmic scale) represent power-law fittings of the relation between L and m. The best-fit exponents in different plots 
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Figure 2. The distributions P(logx(t)) for the particles in our ensembles. Every line of histograms corresponds to a different central 
mass parameter (m = 10“'^, m = 10“®, m = 10“^) whereas every column to a different time snapshot (t = 10®, t = 10^, t = 10®). The 
histograms are shown for e=0.1. 


range in values between p ~ 0.35 and p ~ 0.5. We also note a tendency towards smaller values for smaller t. However, the 
bigger values are more representative of the true exponent, since they appear at times closer to the limit when x(f) tends 
to its limiting value for chaotic orbits, i.e. the Lyapunov charac teristic number (LCN). If I = limf_»ooX(f) denotes the LCN 
limit, it is well known (see e.g. Kalapotharakos fc Voelij 1 20051) 1 that the generic behavior of x(f) is to fall like up to a 
time ti ~ . The time ti is called ‘Lyapunov time’, and represents a saturation time beyond which the curve x(^) starts 

stabilizing towards the limiting value 1. The temporal change of x(t) for times t < ti is reflected in the histograms of Fig[2] 
We can observe that, along a fixed panel row (i.e. fixed m, integration time increasing from left to right), the right wing of 
the histogram is shifted in general to the left as t increases. The shift is more conspicuous in the uppermost panel row (small 
m), in which the orbits have, in general, smaller values of the LCN, and hence, larger values of their saturation times ti. On 
the other hand, in the lowermost panel row (m = 10~^, i.e., large), the saturation time is small (ti < 10® for nearly all orbits). 
Hence, the histogram N{x) remains practically invariant beyond the time t = 10®, as shown in the three panels of the same 


The difference in the saturation times between small and large m has, as a consequence, that in each column of Fig[3](same 
parameters but increasing t), the value of L for all m presents some shift downwards as t increases. The shift is important for 
small values of m, while it is nearly negligible for large values of m (for which the orbits reach their saturation times ti already 
before t = 10®). Hence, the effective logarithmic slope p increases as t increases. Nevertheless, the tendency of p to increase 
with time is only temporary, since p is stabilized after all the orbits have reached their saturation times. This happens around 
t = 10®. 

Another remark is that the dependence of the exponent p on the anharmonicity parameter e appears to be weak. This 
fact confirms that the main source of the chaotic behavior of the orbits is the scattering by the central m ass, while nonlinear 
effects due to the quartic terms in the potential are of small importance. This agrees with findings in i Kandrup fc SiderisI 
1 2 OO 2 I) ). It should be noted, in this respect, that a power-law relation L ~ can be extracted by an a posteriori analysis of 
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Figure 3. The L = versus m relation in logarithmic scale for all our experiments, with parameters as indicated in the panels. 




Figure 4. The L versus m relation as derived from the data of (a) Kandrup and Sideris (2002) (lowest order approximation of a triaxial 
Dehnen potential), and (b) Merritt and Valluri (1996) (potential corresponding to a galaxy with cusp). In (a), the compiled data confirm 
a power law relation L ~ mP over the whole range of values of m considered. The solid line corresponds to a best fit yielding an exponent 
p 0.33. In the right panel (b) the solid and dashed lines show the dependence on m of the logarithms of the mean first and second 
Lyapunov exponents respectively. Both lines extend to a non-zero value of logL for m ^ 0. This represents the values corresponding to 
what is called ‘residual chaos’ in section 4 below, i.e., chaos induced by the central cusp itself. The best-fit exponents were computed in 
this case only by the three rightmost data points. 
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Figure 5. (a) Variations of the value of the formal integral truncated at order 4 (thin line) or 10 (thick line) for an example of box 
orbit with initial conditions x = 0.34, y = 0.23, z = 0.21, Vx = Vy = Vz = 0. (b) The variations of the order 10 truncated series in greater 
detail (notice the change of scale in the vertical axis). 


data i n independent published works i Kandrup fc SiderisI ( 2002 ). Merritt fc Valluri I IQOd )). In particular, Kandrup fc Sideri3 
1 2 OO 2 I) computed Hnite-time Lyapunov Characteristic exponents in a potential representing the lowest expansion terms of a 
Dehnen potential with a superposed softened Keplerian term corresponding to the central mass. From their work (their hgure 
12), the relation be tween L and m can be co mpiled in a log-log scale. As shown in Figl4^, one obtains a power law htting with 
p ~ 0.33. Similarly, Merritt fc Valhiril 1 199(il ) computed the first and second (hnite-time) Lyapunov characteristic numbers of 
chaotic orbits in potentials corresponding to a galaxy with cuspy density proHle and a central core radius. Their results can be 
compared to ours in the limit where the core ra dius (their parameter m p) is larger than the sphere of influence of the central 
mass. This is the case mo = 10“^ in Table 1 of Merritt fc Valluri 1 19961) . compiled in log-log scale in Figlljr. Apart from the 
value m = 0 (corresponding to the horizontal line going to — 00 ), the three available data points appear also to be aligned in 
straight lines indicating power laws both for the hrst and the second Lyapunov exponent. The best-ht exponents are p ~ 0.27 
and p ~ 0.38 respectively, while these values are only indicative due to the scarcity of data and the unknown influence of 
the central cusp of the potential to the result. We note, finally, that the breaking of the power law and the appearance of a 
non-zero value of L as m —^ 0 in FigllJ) is due to a phenomenon called below ‘residual chaos’, i.e. chaos due to the cusp itself. 
This phenomenon is analyzed in section 4. 

As an overall conclusion, a power-law L ~ mF appears quite commonly in numerical computations of the Lyapunov 
exponents of the stellar orbits chaotically scattered by a central mass in various galactic models. We now proceed in a 
theoretical modelling allowing to interpret the origin of this power law. 


3 THEORETICAL MODELLING 
3.1 Transit and out-of-transit dynamics 


As mentioned in the introduction, a modelling of the process of chaotic scattering of the orbits by the central mass becomes 
feasible by considering two distinct regimes of the motion, i.e. i) transits through the sphere of influence and ii) out-of-transit 
oscillations. We begin by showing that within the out-of-transit regime the orbits obey three quasi-integra l s of m otion. Such 
integrals can be written as formal series, using, for example, the ‘third integral’ approach ( ContoDOulos 1 196(]l) ). A formal 
series has the form 4? = >I >2 -|-4?3..., where are polynomial terms of degree r in the canonical variables {x,y,z,px,Py,Pz)- The 
term $2 is set equal to the harmonic energy in one of the three possible degrees of freedom, i.e. {p^ + /2, {py + Uyy'^)/2, 

or (Pz + uJzD)/2. Terms of higher degree are computed recursively, by solving, order by order, the equation 


-^ = [<E>, R 2 + R 4 ] = 0 (10) 

where [•,•] denotes the Poisson bracket operator, and H 2 , H 4 , are the quadratic and quartic terms of the Hamiltonian ([3|). 
This means that the formal integrals are possible to define for the Hamiltonian neglecting the influence of the central mass. 
Then, we test numerically how well they are preserved in the full model, in the out-of-transit regime. 

Equation uni yields, at degree r, the homological equation 
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a) b) 




Figure 6. (a) Time evolution of the quantities $2 for one orbit with energy E = 0.18 in the experiment with m = 10“'^, s = 0.1 

(initial conditions as in FigO. (b) The time evolution of $cc in detail (solid curve, corresponding to the scale of the left vertical axis), 
along with the time evolution of the orbit’s distance from the center (dashed curve, corresponding to the scale of the right vertical axis). 
Gray points on the latter curve correspond to local minima, while gray points on the curve correspond to local extrema. Note that a 
local maximum of the curve ^x{t) coincides always with a local minimum of the curve r{t). 


[H2,^r] = [$r-2,F/4] 


( 11 ) 


We used a computer-algebraic program to solve, step by step, the homological equation, for all three formal integrals defined 
as the above way, up to the 10th degree in the variables {x,y, z,px,Py,Pz)- As an example, for the formal integral tOj, = 
-|- uj'^x'^)/2 + ... up to 4th degree we have 

$ 2 ; = 0.5x^ + 0.5px — £ ^■0025pxPy — O.OOlpxPy — O.OOSAp^pyPz 

— O.OlpxPyPz — 0.00125p^p3 -I- 0.0025pyX^ -|- O.OOSdpyp^x^ 

-I- 0.00125^2®^ — O.OlpxPyxy — 0.007pyXy — 0.0052pxPzxy 

— 0.02pyPxXy - O.OndpxPyV^ + O.OlpxPzy^ - 0.0152®?/® 

— 0.0086pxPyXZ -1- Q-Qlp^xz — 0.005pxPzXZ — 0.0017pxyz 

— O.OApxPyyz -I- 0.0017x^yz — O.OSxy^z — 0.00125p^z^ 

— 0.00125x^2®) -b 0(£^) . (12) 


Similar expressions are found for the formal integrals t&y, $ 2 . The degree of approximation of these expressions can be tested 
by probing how well the integral values are preserved along individual orbits integrated first in the Hamiltonian H = H 2 -b H 4 
(i.e. without the central mass). For orbits in the energy range considered, we find that, up to r = 10, all three integrals 
are preserved to within a time variation of about 10“^^+''/^“^^ at the truncation order r. Figure [5] shows an example of this 
behavior. Panel (a) shows a comparison of the variations of the quantity $ 2 ,, computed as above, at the truncation orders 
r = 4 and r = 10, for an example of box orbit with initial conditions as indicated in the panel. The maximum variation is 
about ±2 X 10““* at the truncation order r = 4, but it reduces to abo ut ±3 x 10~^ at the tr u ncati on order r = 10 (magnified 
in panel (b)). In fact, the estimates of Nekhoroshev theory (see, e.g.. lEfthvmiopoulos et al.l (|200J)) yield that the variations 
continue to decrease up to an optimal truncation order Vopt ~ l/(£i?), in which they become exponentially small, i.e., of order 
0(exp(—l/(£if)). However, even low order truncations are sufficient for estimating the values of the integrals ^x, ^y and <l >2 
in practice. 

Restoring, now, the term —mjlE -bd®)^^® in the potential, we compute the variation of all three quantities ^x,^y,^z 
during both transits and the out-of-transit regime. Figure [6^ shows the example of an orbit with initial conditions as in Fig[5] 
(energy E = 0.18) in the case m = 10”"^, s = 0.1. All three integrals are seen to exhibit abrupt jumps that coincide in time 
(e.g. at the times t = 30, f = 52 or t = 61 in Fig[6K). A closer focus to the jump at f = 52 is shown in (Fig[6]3), superposing 
the variations of $ 2 : in time with the variations of the distance of the orbit from the center. Clearly, the most important jumps 
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Figure 7. Variations of <I>a;(i) for an orbit with energy E = 0.18 and four different values of m (left), or e (right). The values of (/r, e) are 
indicated in the figure. Note the change of scale in the ordinate of the four left panels, indicating that the variations of the quasi-integral 
depend strongly on the value of the mass parameter (while they are nearly insensitive to the non-linearity parameter e as long as the 
latter is not very close to unity). 




Logio crai, LogiQ ni 

Figure 8. (a) The points represent the logarithm of the mean absolute value of the 'hx integral’s jumps logd<I> 2 , for different values of 
the central mass parameter, versus the logarithm of the mean minimum distance from the central mass over all transits for chaotic 

orbits with initial conditions as in Figl6] The fitting straight lines have equation as indicated in the figure, (b) The mean radius rj, within 
which the angular momentum is approximately conserved (see text) versus m. 


take place during transits through the sphere of influence of the central mass (which is about 2 x ~ 0.1 in this case, 

see below, end of subsection 3.1). A similar behavior is found for the quasi-integrals ^y, On the other hand, the values of 
all three integrals are preserved to within two significant figures in the out-of-transit regime (without the central mass, the 
precision increases to about four significant figures at the truncation order r = 4, see above). 

Figure[7l now, compares the variations of ^xit) for an orbit with the same initial conditions but different values of m and 
e. The overall size of the variations appears rather insensitive on e, while the size of all jumps (i.e. the difference between the 
value of fha; at the local maximum and minimum along a jump) clearly increases as m increases. One can observe also that 
the ‘landing’ value of ^x{t) at the end of one jump appears to be more and more unpredictable as m increases. Essentially, 
this uncertainty in the value of ^xit) after each jump is a measure of the chaoticity of the orbit. 
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Figure 9. Upper row: projection of a thin tube orbit around a ‘boxlet’ 3:2 periodic orbit in the planes {x,y) (left panel), and {y,z) 
(middle panel), along with the evolution of its finite-time Lyapunov exponent x(i) (right panel). The central mass parameter is set equal 
to m = 0, and the nonlinearity parameter e = 0.9. The orbit’s initial conditions are: x = 0.1894, y = —0.04796, z = 0.0014, Ua, = —0.3578, 
Vy = —0.4682, Vz = —0.097. The orbital energy is i? = 0.193. Middle and lower rows: same as in the upper row, but with a central mass 
m = 10“'* and m = 10“® respectively. 


Quantifying the behavior of the jumps after many transits allows to see that the scattering of the orbits can be modelled 
essentially by Keplerian hyperbolic dynamics. Figure [8^ shows the mean value of the jumps d'f’a, (measured as the difference 
between the local maximum and minimum values of f&a, at each jump) plotted against the mean value of the radii rmin, where 
rmin (computed from the data of FiglhJj) means the radial distance from the central mass at the point of closest approach 
during one jump. Both means are taken over the ensemble of all jumps during an integration of an orbit up to a time t = 10®. 
The computation is repeated for different central mass values m, keeping the orbit’s initial conditions hxed (same as in Fig|6]). 
During this time, the orbit exhibits some thousands of transits, thus the evaluation of mean values for both and rmin has 
small statistical error. Figure shows the mean value of as a function of the mean value of Tmin, superposing the plots 
in log-log scale for the masses m = 3 x 10“®, 10“*, 3 x 10“*, 10“®, and 3 x 10“®. The straight fitting lines have inclination 
-1, whereas we note that the vertical distance between two successive lines is about equal to log^Q 3 « 0.5. Thus, the mean 
jump d^x is consistent with a (absolute) Keplerian potential law: 

= . (13) 

min 

The above result concerns the variations of the quasi-integrals valid in the out-of-transit regime. On the other hand, inside 
the sphere of influence the force field is approximately central, thus another approximate integral, of the angular momentum 
L, is valid. During transits, we follow the time variations of L and determine, at each transit separately, the maximum radius 
tl up to which the variations AL are smaller in size than a fixed percentage of L (taken as 10% difference measured from the 
value of L at the point of closest approach to the central mass). Figure[8]D shows the so-obtained mean value of vl as function 
of the central mass parameter m in log-log scale for the same orbits as in FiglH^. Albeit with considerable scatter, the data 
allow to determine a best-fit power law ~ 2.2m°'®^®. This is close to the power law for the sphere of influence (Eq.®, 
thus allowing to identify tl as a measure of Vm yielding the estimate rc ~ 2. We note that the threshold of 10% variation 
of the angular momentum is rather arbitrary. However, it allows to obtain a meaningful result for box orbits which, far from 
the BH sphere of influence, undergo angular momentum variations of order 100% (with more stringent thresholds we identify 
significantly fewer transitions through the BH sphere of influence). 

In conclusion, their scatter notwithstanding, the previous plots suggest a qualitative picture in which the dynamics of 
chaotically scattered centrophylic orbits can be modelled as a sequence of i) transits through the sphere of influence, in 
which the orbits follow approximately a Keplerian hyperbolic dynamics, followed by ii) box-like wanderings within the rest 
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Figure 10. (a) Comparison of the time evolution of the stretching number a{t) (Eq. ||14I| ~) with the variations of the quasi-integrals $ 3 ;, 
^y, (b) The mean stretching number a after thousands of transits versus m for the orbit with same initial conditions as in Fig[ 6 l A 

power-law fit yields a oc mP'^^. The theoretical estimate a oc is found in the appendix. 


of the available space in the interior of the equipotential ellipsoid corresponding to a hxed value of the orbital energy. This 
picture is quite generic when the three frequencies wi, LJ 2 , ws are far from low-order commensurabilities. The generation of 
such commensurabilities necessitates a se parate treatment, since, then, the orbital sample contains also many resonant boxlets 
I Miralda-Escude fc Schwarzschild (198^)). In onr model, boxlets corresponding to low-order commensurabilities are generated 
by the qnartic potential terms in Eq. for large (0(1)) values of e. An example is given in Fig[9l for e = 0.9. When m = 0 
(top row) the orbit with initial conditions as indicated in the figure’s caption is a three-dimensional thin-t nbe orbit around a 
resona nt 3 : 2 orbit which exists in the plane {x,y). The periodic orbit is stable, and, hence, centrophobic lIMerritt fc Valluri 
1 19^ 1. In FigH top row, the tube orbit around the boxlet also avoids the center. Hence, it is a regular orbit, as confirmed 
by computing its Lyapunov characteristic exponent (right panel), which goes to zero. Now, the orbit’s closest approach to 
the center is at a distance r ~ 2 x 10~^. Thus, by adding a central mass with parameter m = 10““* (middle row), the orbit 
now crosses the central mass sphere of influence (of radius ~ 10~^^®). As a result, we observe that the orbit is significantly 
deformed, and looses its resonant character, while the Lyapunov exponent stabilizes to a positive value ~ 10~^, i.e., the orbit 
becomes weakly chaotic. For still larger m (m = 10~^), the orbit exhibits the usual behavior of a chaotic centrophylic orbit, 
with a Lyapunov exponent ~ 2.5 x 10~^. 

We now model the chaotic scattering process of centrophylic orbits in order to derive theoretical estimates for the orbits’ 
Lyapunov exponents. 


3.2 Theoretical Estimates on Lyapunov Exponents 


Consider an orbit integrated along with the variational vector ^{t) of a nearby orbit np to the time t. Let be the modulus 
of the deviation vector of the orbit at the ti me tj = iAt, where At — t/n is the timestep corresponding to a splitting of the 
integration in n steps. The stretching number Voglis fc ContopoulosI (1994) at the i-th step is defined as 


1 , 

The finite time Lyapunov characteristic number is equal to the mean stretching number along the orbit, since 

n 

1 V- 


t ^0 nAt -4^ Ci-i 

i = l 


— y ai 

i=l 


(14) 


(15) 


Figure [TOk shows the time evolution of the stretching number a(t) as a function of time, along with the variations of the 
integral ^x{t) for the orbit with initial conditions as in Fig[71 for m = 10 ~®, e = 0 . 1 . Far from transits, the function a{t) 
shows an oscillatory behavior around zero. This behavior is characteristic of a nearly-harmonic oscillation, while the quartic 
potential term implies an overall linear growth of the deviation vector in the out-of-transit regime, which is of order 0{e)t. 
On the other hand, in Fig IlOb the curve a{t) clearly exhibits positive peaks at all transits. 

Consider the set S = {ii, * 2 , inr} such that the orbit is in transit at the time ti with i £ S, nr denoting the total 
number of timesteps during which the orbit is in the transit phase. Let S be the complement of S with respect to the set 
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{1, 2,n}. One has the estimate: 


nAt ^ O.-i 


o 


( In(enAt) 
I nAt 


ij — l 

implying that the contribution of the ‘out-of-transit’ stretching numbers to the final value of x(i) goes to zero as ~ Int/t. 
Thus, one has the approximation 

ri'j' 

ij GSJ=1 

Setting a = (I/tit) we get x(^) ~ {nT/nAt)a. The quantity 

N ■ = dLE 

""" nAt 

hereafter called ‘rate or visits’, represents the number of transits per unit period of an orbit whithin the sphere of influence. 
We then have 


X{t) « Nvisa 


(16) 


Estimates on the mean value of x(^) for all transiting orbits at fixed central mass parameter m will then follow by estimating 
separately the quantities a and Wis- 

Assuming, as evidenced above, that the transits are governed by nearly-Keplerian hyperbolic dynamics, in the Appendix 
it is shown that for orbits of given energy E, one has the theoretical estimate 

• (17) 


Figure fTUb shows the mean a computed numerically for an orbit of fixed energy with the same initial conditions as in Fig|^ 
integrated under various values of m. Numerically we find the exponent 0.28, which is in fair agreement with the theoretical 
exponent 1/3 of Eg. (1171) . The predicted dependence of o on the energy, probed numerically below, is also verified. 

We now focus on estimating N^is. The frequency whereby an orbit visits the sphere Vm depends on the geometry of the 
orbit in configuration space. We distinguish two cases, explained with the help of Fig llll 

i) 3D-orbit: as long as the three quantities $ 2 ,, $2 obtain comparable values, an orbit fills nearly uniformly the 

available conhguration space, which has the form of a deformed 3D box. 

ii) planar orbit at least one of the three quantities $ 2 ,, ^y, > 1>2 obtains a value smaller than a given threshold (given 
by Eg. (1181) below). Geometrically, the amplitude of oscillations in at least one of the three axes in the out-of-transit regime 
should be smaller than the radius r^n of the sphere of influence IFig lllh . schematic). Quantitatively 

< rm (18) 


where k stands for x, y, or 2 . 

Note that ‘linear’ orbits, i.e., tubes around the stable axial orbits also exist, but their importance is rather limited because 
they are considerably fewer than the planar or 3D orbits. 

The dependence of the rate of visits to the central masses’ sphere of influence on the geometry of orbits can be modelled in 
the following way: for 3D-orbits, considering all possible straight line segments connecting two different points on the surface 
of the box delimiting the orbit (see Fig lllh i. N^ia can be approximated as proportional to the percentage or line segments 
passing through the sphere of influence. Then, N^ia oc Sr^/Stot, where Sr^ S'UcI Stot are the surface of the sphere of influence 
and of the box respectively. The linear dimension of the box is of order I ~ E^^^, where E is the energy of the orbit. Thus, 
Wis oc rf,/E, or (taking into account Eq.((5|) 

rn^!^ 

Wis oc . (19) 


For planar orbits, one has, instead, the estimate oc r^/i, or 


N' ■ 


2I/3 


£ 1/2 £( 1/2 


( 20 ) 


These estimates are confirmed numerically. Fig lllb shows a computation of the rate of visits Wis, A/is for two orbits with 
constant energy E — 0.18, but for different mass parameters m. The number of visits within a total integration time t = 10® 
are counted, and criterion dH} is used in order to distinguish either orbit as ‘planar’ or 3D (the corresponding rate of visits is 
found by dividing the number of visits by t = 10®). The difference in the shape of the orbits is evident fFies Illb -f'). Returning 
to Fig lllb . the best-fit exponents of the relations Nvia and NLa to m are 0.7 and 0.34 respectively, in fair agreement with 
Fas. (1191) and (1201) for constant E. 

Taking into account eqs. (I20J,(IT91) and 113 we find for the Lyapunov number of 3D-orbits of energy E the estimate: 
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Figure 11. (a) Schematic representation of the configuration space of motion of “3D-orbits” and “planar” orbits. The cube represents a 
box volume covered by a “3D-orbit”. The central sphere represents the central mass sphere of influence. The three thin parallelepipeds 
(gray) correspond to one parallelepiped side being equal to 2rm- Planar orbits are orbits lying inside one of the three parallelepipeds, 
(b) The number of visits (up to time t = 10^) versus m for two orbits characterized as ‘3D’ (black points, initial conditions x = 0.227, 
y = 0.155, 2 : = 0.139, Vx = Vy = Vz = 0), and ‘planar’ (gray points, initial conditions x = 0.227, y = 0.155, 2 = 0.0,Ucc = Vy = Vz = 0), 
for e = 0.1. The straight lines represent power-law fittings yielding the best-fit exponents 0.7 and 0.34 respectively, (c-d) Projection of 
the ‘3D’ orbit in the planes {x,y) and (x, z). (e-f) Same for the ‘planar’ orbit. Note in (e) the change in morphology induced by a big 
jump in the values of the quasi-integrals Such jumps are stochastic in nature, and they may occasionally convert a 3D orbit to 

planar and vice versa. 


E ^ £ 1/2 ~ £ 3/2 
whereas, for planar orbits 

, ^ rrd!^ 

X ~ ^1/2 ^ ^1/2 ~ E 


( 21 ) 


( 22 ) 


Figure [12^ shows the values of logx against log£ for all the orbits in our considered ensembles for six different values of 
the mass parameter m as indicated in the caption. The varions ensembles are clearly distingnished by the concentration of 
points in the scatter plot, the uppermost concentration corresponding to the ensemble in the experiment with the highest 
mass parameter (m = 10“^). Each ensemble can be roughly described as consisting of a ‘rising’ and a ‘falling’ part of the 
valne of logx vs. log£. The two parts meet at a point of maximnm of logx- The position of the maximum moves to the 
right with respect to log£ as m increases. However, the level value of logx fhe maximum remains remarkably constant, 
i.e. nearly independent of m. 

The straight lines show power law fittings of x with £ for the falling part. Despite the large scatter of the data points, 
we find indicative logarithmic slopes lying in the range between -1 and -1.5 for all ensembles considered. This behavior will 
be explained below. A power-law roughly holds also in the rising part. The point of maximum corresponds to abont the point 
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Log,o E 

Figure 12. (a) The logarithms of finite-time Lyapunov characteristic numbers logx of all orbits in the choosen ensembles versus the 
logarithm of the orbital energy log£^ for the experiments with mass parameters (from top to bottom) m = 10“^, m = 3 x 
m = 10“^, m = 3 X m = 10““^ and m = 3 x 10“^, at the end of the integration {t = 10^). The straight lines represent asymptotic 

power-law fittings for the right parts of the plots separately for each mass parameter, (b) The energy Emax where logx in (a) exhibits 
a global maximum versus m. The power-law fiting is Emax = (c) The plot logx versus = logE in greater detail for the masses 

m = 3 X 10~^, and m = 3 x 10”'^. In each case, numerical values are distributed between two lines with inclination -1 and -1.5, as 
predicted from Eqs. 112211 and lETJ for the planar and the 3D orbits respectively. 


where the associated best-fit power laws intersect. The intersection point defines an energy E^. Computing and plotting Ec 
against m yields approximately a power-law Ec oc ('Fig |12b 'l. 

These features can be understood by the following consider ations: first, one can note that the left part represents regular 
or sticky chaotic orbits which have the morphology of pyramids Jlvlerritt fc Vasilis (2011)), i.e. they lie nearly entirely within 
the sphere of influence of the central mass. Such orbits can be described by perturbations to the Keplerian dynamics of the 
central mass. The limiting energy value Ei up to which an orbit lies entirely within the radius r = Vm can be estimated by 
requiring that the sphere r = Vm constitutes the surface of zero-velocity. The estimate 


A V m A, 


.2/3 


1 


-) m 

r, 


2/3 


(23) 


holds, where fl is a geometric-mean estimate of the harmonic frequencies, bounded by the highest of the three frequencies 
u)x,^y,u)z. The dependence of Ei on m in Eg. (1231) has the exponent 0.66, close to the exponent in the numerical htting of Ec 
versus m ('Fig ll2b 'l. This suggests that Ei ~ Ec (the near equality holds also checking numerical vs. theoretical coefficients). 
Use of the estimate (1231) is made in the next subsection. 

On the other hand, most chaotic orbits lie beyond the energy Emax- In fact, isolating the plots logx logiS (Fig ll2b 'l 
for two values of the mass parameter m allows to see that the whole ensemble of orbits in the right wing are delimited between 
two limiting lines with inclinations -1 and -1.5 respectively, i.e. as predicted from the estimates of Eqs. (1221) and (EH for the 
planar and the 3D orbits respectively. The coexistence of ‘planar’ and ‘3D’ orbits explains in this way the scatter in the data 
points of Fig |12h . 
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Figure 13. (a)Mean finite Lyapunov characteristic number L vs. the mass parameter m as it follows theoretically from Ea. (l28l l (black 
dots), compared to the numerical computation (gray dots) at the integration time t = 10^, for all the experiments with e = 0.1. The 
mean inclination is close to 0.5 in the range 10~^ < m < 10~^. (b) Theoretical computation of L versus m in a ‘isothermal’ number 

1 /2 

density model (see text). The triangles and gray fitting line correspond to the choice tr = Vq ' /4, filled circles and black solid line to 
cr = V’jj^^^/2, squares and dashed line to tr = . 


3.3 Final theoretical estimates: the power law L oc mP 


Assuming (as evidenced in FigITJo) that a time t = 10® is sufficient for a saturation of x(f) close to the limiting value of 
chaotic orbits, i.e. close to the Luapunov characteristic number (LCN), the mean LCN of the orbits in an energy range 
Emin Si F ^ Emax Can be estimated as 


1 

'Tfo 


N{E)x{t = 10®,F)dF 


(24) 


where No = N{E)dE, N{E) is the number density of orbits of energy E, and y(t = 10®, E) is a mean estimate of the 

^min 

value of X for orbits of energy E at the final integration time. Considering only transiting orbits, we set Emin = Ei ~ 
(Ea. (l23ll ). Also, from Fig |12h it is clear that the maximum energy Emax = 0.2 considered in our samples is sufficiently high for 
the orbits’ values of x to fall with respect to the maximum by at least one order of magnitude in the worst case (m = 10“^), 
and typically by several orders of magnitude. 

The mean value x{t = 10®, E) can be now estimated by considering the separate mean values of x S'® well as an E- 
dependent varying proportion of planar vs. 3D orbits. The percentage A of the planar orbits is estimated by the ratio of the 
surface occupied by initial conditions of planar orbits on the box surface corresponding to an energy level E (see Fig lllh ) 
over the total area of the box. Thus 
2ml/® 


A = 


S2d 


Stot 24F Fi/2 • 

The percentage of 3D orbits is 1 — A ss 1 — 2mi/®/Fi/^. Finally, X is estimated according to Eqs. dm), 

_ _ m 

Xp/anar(-^) ^ ^ ’ X3 d(-^) ^ 

with Cl and C 2 being constants of order unity. Then, for the mean value of x s,t fixed energy we have the estimate 


(25) 


(26) 


X{E) « ^XplanariE) + (1 “ DX3d{E) 


(27) 


The main uncertainty in Ea. (l24ll regards the form of the number density function N{E). In self-consistent models, N{E) 
is determined by the distribution function of the centrophylic orbits. On the contrary, in ‘ad hoc’ potential models, N{E) 
cannot be determined self-consistently unless one possesses information on the kinematic distributions allowing to solve the 
reverse problem density —^ N{E). Assuming no detailed model, we hereby estimate the integral of Eg. (1241 using two different 
estimates of N{E) as follows: 

i) we consider the case of a nearly uniform distribution N{E) = const. Combining Eos. (1241 . (1271 and (1261 we obtain: 


L 


0.2 — 


[4cim^/® -I- 10c2m'i/® 


2.23(4ci -I- 2c2)m] 


(28) 
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The exact dependence of L on m in the model (I28II depends on the relative values of the constants ci and C 2 , as well as 
on factors entering in all the above estimates. Figure [13^ shows L against m in logarithmic scale, by the estimate (I28II (black 
points) setting simply ci = C 2 = 1. The plot indicates an approximate power-law L ~ mX. This is produced as follows: since m 
is a small quantity, the leading term in Ea. (l28ll is the one with lowest exponent, i.e., . Thus, in the absence of additional 

terms, we would have p = 2/3. However, for relatively large m, the second most important term (linear in m) has a negative 
sign. Thus, it lowers the rightmost part of the curve L vs. m (the presence of in the denominator only marginally affects 
the overall power-law behavior). If we approximate the new curve by a single power-law fitting, we then find a lowering of the 
exponent p, i.e. 

L oc (29) 


with q varying between 0.1 and 0.2. However, one notices that the whole curve in log-log scale deviates downwards from the 
power-law approximation (|29l) for the highest mass parameter values, i.e. q has a weak (increasing) dependence on increasing 
m. Comparison with the numerical data (gray dots, for e = 0.1) shows that the model reproduces the slope of the numerical 
curve, which also exhibits a lowering of the value of L at high mass parameter values with respect to an exact power law. 
Overall, the theoretical curve has a factor ~ 2 difference from the numerical curve, which is consistent with uncertainties in 
the theoretical coefficients. Notice also that a systematic lowering of the values of L with respect to a power-law fitting is 
discernible in all the plots of FigO Finally, we have checked that the appearance of an approximate power-law persists, with 
exponents around p « 0.5, for different choices of the constants Ci and C 2 . This shows that the dependence of the integral 
on m (which enters in the integral as a parameter) is not sensitive on the details of the distribution of the planar vs. 3D 

orbits. _ 

ii) Fig |13b shows an evaluation of the integral (I24II for an isothermal (or ‘ergodic’, see Binnev fc Tremaine! ( 2008h l model 
of N{E), i.e. N{E) oc . The constant cr is a measure of the velocity dispersion in the central parts of the galaxy. 

Assuming a core density p « 3/(47r) (in our units, corresponding to total mass M = 1 at radius R — 1), by the Virial theorem 
has to be taken of the order of the absolute value of the central potential well Vo ~ VixGprdr ~ 1.5. Figure [TSjr shows 
the evaluation of the integral (1241) for three different choices of cr, namely cr = a — (1/2)Iq^^^ and cr = (l/d)!/,^^^. In all 

three cases we recover here as well an effective power-law behavior, with not very different exponents, i.e. p ~ 0.58, p ~ 0.56 
and p ~ 0.56 respectively. 

In conclusion, we find that an approximate power-law relation of the form (12911 is robust against details of the form of 
the function N{E). 


4 EFFECT OF CENTRAL CUSP 


In the model ®, the existence of many box orbits was a priori guaranteed due to the harmonic core in the center. It is well 
known, h owever, that realistic mode ls o f the cen t ral pa rts of galaxies include central density cusps p{r) ~ r ' (see e.g. the 


review m 


Binne y fc Merrifieldl ( 199^1. or iMerrittI (IiQQqIH. In su ch models, the cusp itself transforms most centrophylic orbits 


to chaotic ( Merritt fc Valluri ( 1996ll . Merritt fc Quinlan ( IQQsI H. Even without central black hole, one then expects the orbits 


to exhibit positive Lyapunov exponents. We hereafter call this effect ‘residual chaos’, i.e. chaos existing even when m = 0. 
The corresponding mean Lyapunov exponent of the centrophylic orbits is denoted by Lq. 

Adding, now, a central mass (m yf 0) we seek to determine the dependence of L on m. The theoretical analysis of the 
previous section formally breaks down, since one cannot define the formal integrals 4>3;, 4?^, <I >2 even for m = 0. We thus rely 
on numerical computations. To this end, we consider again the Hamiltonian function ©, changing the potential model to 

7T) 

( 30 ) 

where Vd represents the ellipsoidal Dehnen model ( Dehnen 1 1993li i: 

[V’(oo) — i/’('w)]dr 


VD{x,y,z) = —TvGabc 


where 


•^(t -I- a^)(r + iP){T + c^) 


(31) 


with 


and p{w) given by 


ij){w) = / p{w'^)dw'' 

Jo 


w = ^ + a ^ b ^ c > 0 

a‘‘ o 2 


p(y;) = - 0)^ W ^(l-fui) , 0^7 <3 . 

47raoc 


( 32 ) 
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Figure 14. The mean Lyapunov exponent versus the central mass parameter m in four Dehnen triaxial models with central cusps (see 
text). Black filled circles correspond to the data points for the cusp exponent 7 = 0.3, gray triangles to 7 = 0.7, black rombuses to 
7 = 1.3, and gray squares to 7 = 1.7. The straight lines (black or gray) are power law fittings obtained by the rightmost seven points 
for 7 = 0.3, five points for 7 = 0.7, five points for 7 = 1.3, and four points for 7 = 1.7. The horizontal dashed lines correspond the the 
values logj^pLo, where Lq is the ‘residual Lyapunov exponent’ in each case (see text). The points where each horizontal line intersects 
the fitting line for the same corresponding 7 are marked as stars. 


The parameters a, b and c correspond to the lengths of the major, intermediate and minor axis of the triaxial equipotential 
surface corresponding to w = 1. The parameter M determines the system’s total mass. We use a similar algorithm as in 
Merritt fc Fridman ( 1996h in order to numerically evaluate the integral (I31II as well as its spatial derivatives, i.e., the forces. 


A ‘weak cusp’ corresponds to 7 < 1. In this case, the modulus of any component of the force generated by Vd goes to 
zero at the center, reaches a maximum value at a certain distance from the center, and then falls-off tending to zero at large 
distances by a Keplerian law. This behavior of the force allows to determine values of the parameters a, b, c, and M, for fixed 
7 , so as to create a system exhibiting a similar geometry and value of the total mass as the simplified system corresponding 
to the potential O of section 2 , the two systems being, hence, differentiated essentially only by the presence of a central cusp 
as opposed to a harmonic core respectively. The parameter determination is realized by the following algorithm: 

i) We set the ratios a : b = ujy : oJx and a ■. c = uJz/ijJx, where ujx, ojy, and lOz are the parameters of Q. 

ii) We fix the value of a so that the quantity dVoldx presents maximum at the point {x = Xmax, y — 0, z = 0), with Xmax 
choosen so as to represent the point where the harmonic model in ([4| yields a total mass equal to unity. We find Xmax — 1.06. 

iii) We fix M so that the force Fx under the potential Vd be equal to the force Fx, under the potential with m = 0 , at 
the point (x = Xmax,y = 0, 2 = 0). Note that, since the value of the force depends essentially only on the total mass inside a 
given radius, this normalization means also that the total mass inside an ellipsoidal surface crossing x = Xmax is nearly equal 
in the harmonic and in the 7 -models. 

We examine two values of 7 in the weak cusp case, namely 7 = 0.3 and 7 = 0.7. 

In the case, now, of a ‘strong cusp’ (7 > 1), criterion (ii) can no longer be implemented, since dVoldx —>■ 00 as a; ^ 0 
(with j/ = 2 = 0 ), implying that the quantity dVoldx does not present any smooth maximum along any of the principal 
axes. As a simple (but somewhat arbitrary) way to bypass this difficulty, we keep a constant to the value a = 6.67 found by 
criterion (ii) in the second ‘weak cusp’ experiment (7 = 0.7). Then, we fix the remaining constants by criteria (i) and (iii). 
We run also two strong cusp experiments, with 7 = 1.3 and 7 = 1.7. 

In all four experiments, the initial conditions are choosen as in section 2, namely 200 initial points of zero velocity 
randomly distributed on equipotential surfaces with V = E and E choosen uniformly in the range 0 ^ ^ Emax, with Emax 

choosen as Emax = Voixmax), so as to ensure that the resulting centrophylic orbits exhibit oscillations of amplitude at most 
equal to Xmax- However, here we integrate the orbits only up to t = 1000, since the complexity of force evaluation in the 
model Vd renders the computational cost of longer integration prohibitive. Yet, as shown below, our smallest found Lyapunov 
exponents are about L ~ 10~^ ®, implying that a time t = 10® is marginally greater than the saturation time t ~ XjL even 
for the orbits with smallest Lyapunov exponents. 

Figure [ 14 ] (analogous to FigO top row) shows the mean Lyapunov number L — x{t = 10®) for our ensembles of orbits 
in the four above experiments, as a function of the central mass parameter m. We note immediately that power-law fittings 
are possible in only a range of values of m, i.e. above a critical threshold value m > mc( 7 ). In Fig |14l an estimate of the 
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threshold value is found by the abcissa of the points displayed by stars. They are computed as follows: the four inclined lines 
represent power-law fittings for the rightmost part of the numerical curve of L vs. m in each experiment. The horizontal lines 
illustrate the level values of the quantity Lq = L{rn = 0). We call Lq the ‘residual Lyapunov exponent’. It represents the mean 
Lyapunov exponent of the centrophylic orbits when m = 0, i.e., under the influence of the central cusp only. The point at 
which a horizontal line of fixed Lo intersects the corresponding inclined fitting line of L vs. m in the same experiment (same 7 ) 
marks the position of a star-point, and the associated abscissa, i.e. a critical mass value me From Fig ll4l it is straightforward 
to see that both Lq and rric increase in general with the strength of the cusp (i.e. the value of 7 ). On the other hand, it is clear 
from the numerical data points that an approximate power-law correlation between L and m persists, in all four experiments, 
for central mass parameters larger than m = m^. A physical understanding of this phenomenon is the following: the chaotic 
scattering caused by the cusp itself acts dynamically as a central mass concentration, whose distribution is not point-like but 
follows the cusp density law. As long as the BH mass is small, the effect of the cusp is dominant over the efect of the BH. 
Thus, the mean Lyapunov exponent of the centrophilic orbits remains nearly equal to the residual Lyapunov exponent Lq. 
But beyond the BH mass scale m = rric, the black hole dominates over the central cusp. Then, we recover a correlation of L 
with m. This goes asymptotically to an effective power-law. Furthermore, the exponents found by fitting in the range m > rric 
are all about p ~ 0.4, i.e., not very different from those of the corresponding data in FigO the essential difference in the two 
plots being with respect to whether or not we observe a a critical mass scale in which the power-law br eaks. 

We note finally, that despite the sparsity of their datapoints, the results of Merritt fc Valluri (199^), reproduced here as 
FigHbi show essentially the same st ructure as the results of Fig ll4l Thus, the residual chaos phenomenon explains the plateau 
of the curve L vs. m in the data of Merritt fc Valluri (199^) as well. 


5 THE LexM^ LAW IN DISC-BARRED GALAXIES 

N-body simulations of barred galaxies (e.g. Friedli fc Benz ( 199.iI L Friedlil 1 1994 L Norman et al. (199^)) have demonstrated 
that the growth of a central mass concentration induces secular evolution in such systems as well. In fact, although dy¬ 
namically not favored, black hole growth to a mass level as high as IO^Mq-IO^Mq, corresponding to a few percent of the 
mass of a typical galactic bar, could induce even a total destruction of th e bar, with its c o nversion into a nearly axisym - 
metric bulge - like c o mponent. Test particle i ntegrations in bar r ed po tentials (Pfennigei ( 1984 L Pfenniger fc de Zeeuw (198^), 
Hasan et al. ( 199,ll 'lL Norman et al. (199^), Shen fc Sellwooci I 2004I L indicate that a primary mechanism responsible for the 


secular evolution of bars, and even bar dissolution, is chaos induced by the central mass. 

Hereafter we study the dependence of Lyapunov exponents on the central mass paramerer in rotating disc-barred galaxies. 
Two points should be immediately emphasized: i) our modeling in previous sections was based on the existence of box¬ 
like centrophilic chaotic orbits. Such orbits cannot exist in rotating disc-barred galaxies. However, as shown in the sequel, 
centrophilic orbits appear around the main families of planar periodic orbits (e.g. the xi family). Note that the presence of 
some type of centrophilic orbits is an indispensible feature of bars with a rising density profile in the center. As discussed below, 
albeit different in morphology, the centrophilic orbits in barred galaxies are found numerically to exhibit a similar chaotic 
behavior as the boxy centrophilic orbits in elliptical galaxies , ii) Besides the central mass, chaos is generated by the interaction 
of resonances in the corotation domain JContopouloa ( 19^ L Pfennigei lll984ll. S2Mke_&SelIwood[j[l987 Pfenmgei_&FHe^ 


1 199lh . Kaufmann fc Contopoulos ( 1996 I l Patsis et al. ( 1997 1. FuxI ( 200lh . lPichardo et aL ( 2004 1. Kaufmann fc Patsij 1 200!Th l. 
Nevertheless, this type of chaos is a quite distinct phenome non. In fact, most chaotic orbits in the corotation domain belong 
to the so-called ‘hot population’ ( Soarke fc Sellwood ( IOStI II. hence they are not centrophilic. 


5.1 Potential model 

As a case study, we employ the barred-galaxy potential introduced by Kaufmann fc Contopoul^ ( 199^ 1 in a rough self- 
consistent modelling of the galaxy NGC3992. Adding a component for the central mass, the total potential is analyzed as: 

Vtot = Vbh + Vli + Vd + 14 (33) 

where Vbh is the potential generated from the central mass (black hole), while 14, 14, 14 are dark halo, disc, and bar potential 
components respectively. The potential of the central mass is, as before, 
rribh 


\/r'^ + d? 

The remaining terms are as in Kaufmann fc Contopoul^ (199^). The dark halo term is 
I4(r) = , . 

The disc term corresponds to an exponential disc: 


(34) 


(35) 
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Bar : 

Mb 

Q 

b c Dp 


1.5 

5.5 

2.1 0.55 43.6 

Disc : 

So 

Ed 



750 

0.235 


Halo : 

Mh 

bh 



27.5 

12 



Table 1. Parameters of the disc-barred galactic potential. The units are Kpc ^ for Sd-, Kpc for a, b, c and bh, for Mh-, Mh, 

Kms~^Kpc~^ for fip, M^jpc? for Eq. 


V^ir) = -So^r[7o(le,r)i^i(l6,r) - 7i(ie,r)7^o(le,r)] (36) 

where 7o, 7i and Kq, K\ are modified bessel functions of the first and second kind respectively. The bar term is of the Ferrers 
n = 2 type, with the major axis alligned with the y-axis: 

Vdx. y, z) = -^^^[3{2Wiiox^y^ - Wi 2 ox^y^ - W 2 iox^y^ - Wwoy^ + Wo 2 ox^ + 

96 

-VFoioa:^) + VFooo - VFoso*® - VFaooJ/®] (37) 


where th e coefficients Wnk are given by ell iptic integrals. All model’s parameters, as well as the value of the pattern speed Q.p 


Kaufmann fc Contopoulos ( 1996l l. referring to the model for the galaxy NGC3992. They are summarized in table I 


below. Note that the original model contains also a spiral-arm term, which, however, is only important at radii beyond the 
end of the bar, and it is here ignored. 


5.2 Numerical experiments 

As in section 2, we numerically integrate the equations of motion, as well as the variational equations, for planar orbits under 
the Hamiltonian (in cylindrical coordinates): 

H{r,e,pr,pe) = i(Pr + ^) - %Pe + V{r,e) = Ej (38) 

where Ej is the Jacobi constant. The Hamiltonian (1381) describes the motion in a rotating frame with pattern speed flp, while 
pe is the angular momentum in the inertial frame of reference. 

Initial conditions are choosen in a way so as to ensure that they give rise to centrophilic orbits. To this end, we consider, 
as above, ensembles of 200 orbits with initial positions uniformly distributed on a cycle of radius r = 0.1 around the galactic 
center. Initial velocities are given in the direction radially outwards, with modulus choosen so that the value of the Jacobi 
constant is uniformly distributed in the range —2.16 x 10® ^ Ej ^ —2.03 x 10®. This range is choosen so as to correspond 
to energies well below the value at the Lagrangian equilibrium point Li, i.e. 75j,i = —1.915 x 10®. The corresponding orbits 
lie then always inside the corrotation domain, i.e. they support the bar. Orbit ensemble integrations are done for a time 
t = 10® (in comparison, orbital periods are of order ~ 0.1). In different experiments the central mass varies in the range 
10 ® Mq ^ rrihh ^ 10 '^ Mq. 

Figure [15] shows the main result: the mean Lyapunov number L of the chaotic orbits choosen as above is plotted against 
the mass parameter m, choosen here as the ratio m = mhhlrnbar, since this ratio is relevant to a quantification of the rate of 
secular evolution of the bar. We observe again that the numerical data can be fitted by a power-law L oc mP, with p = 0.51. 
We now interpret the mechanisms of chaos and identify the families of orbits which are responsible for this behavior. 


5.3 Interpretation 

The mechanism by which a central black hole generates chaos in a disc-barred galaxy can be visualized with the help of phase 
portraits, obtained by means of a suitable surface of section. Here we employ the apocentric condition r = 0, r < 0, in order 
to define the surface of section. We then plot the intersection points of all orbits with the above section as projected on the 
plane {0,pe). 

Figure [16] shows the surface of section portrait at the energies (Jacobi constant values) Ej = —204000 and Ej = —195000, 
without central black hole lFig ll6h .dL or with a black hole of mass rubh = 10®M© fFig llOb .el. The change of phase space 
structure is evident, namely the insertion of the central mass destroys many rotational KAM curves, corresponding to regular 
(quasi-periodic) orbits around the galactic center. Also, at the second energy level (lower panels), which corresponds to motion 
closer to corotation, a number of librational KAM curves around the 1:4 island of stability are destroyed. In fact, one can 
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Logio ai 


Figure 15. Logarithm of the mean Lyapunov number of the orbits in the ensembles of initial conditions as described in the text versus 
the central mass parameter m = (in logarithmic scale) in a disc-barred galaxy model. The logarithmic slope indicates a 

power-law exponent p = 0.51. 


see that at both energy levels the orbits converted to chaotic are tube orbits around the 2:1 and 4:1 branches of the xi stable 
periodic orbit. The periodic orbits themselves at the corresponding energies are shown in Fig |16b and f respectively. The 2:1 
orbits exist up to energies « —199000. When mbh = 0, the quasi-periodic tube orbits around a 2:1 orbit give rise to two islands 
of stability in the surface of section. The outermost librational invariant curves of these islands correspond to thick tube orbits 
fFig |17l) . The crucial difference between the two tube orbits in Figs ll7b , and b is that in Fig |17b , the tube orbit leaves a hole in 
the center, whose dimension is larger than the black hole’s sphere of influence. On the contrary, the tube thickness of the orbit 
in Fig |17b is larger than the half-width of the orbit’s amplitude of oscilation in the y-axis, i.e., the orbits leaves no hole in the 
center. Then, after the insertion of the central mass, the orbit with same initial conditions as in Fig |17b . retains its regular 
(quasi-periodic) character (Fig |17b ). while the orbit with same initial conditions as in Fig |17b becomes chaotic (Fig |17b l). A 
similar criterion applies to whether or not a thick tube orbit around the 4:1 periodic orbit becomes regular or chaotic after 
the insertion of the central mass (Fig |17b -h). In fact, one readily finds that the initial conditions separating these two types 
of orbits correspond to the last librational KAM curve in the 4:1 island of stability of the section of Fig llGti . 

Investigating the efficiency of the above mechanism to other resonances, one Hnds that the mechanism is not able to 
produce chaotic orbits in the cases of other low-order resonances like 3 : 1 and 6 : 1. In fact, the tube orbits trapped around 
these central periodic orbits in this model form rather small islands of stability. Thus, the tube thickness is small, and we hnd 
no tube orbits able to cross the central masses’ sphere of influence. A similar restriction holds for higher order families in the 
same model. 


6 CONCLUSIONS 


In the present paper we analyze the origin of a numerically observed approximate power-law relation L oc m^, with p ~0.3-0.5, 
where L is the mean Lyapunov exponent of centrophilic orbits in galaxies with central masses (black holes), and m the mass 
parameter, i.e., the ratio of the central mass over the mass of the galaxy. Also, we find that such a law can be recovered in 
quite different contexts and models of galactic systems, ranging from elliptical galaxies with cores or cusps to rotating barred 
galaxies. In particular: 

i) We first make numerical experiments with a simple model of elliptical galaxy with smooth central force field, to which 
the force held of the central mass is superposed. The experiments conhrm the power-law L oc mP, when L is estimated through 
its ‘hnite time’ analog, i.e. the mean value of hnite-time Lyapunov exponents. We demonstrate the statistics of these values 
for centrophilic orbits. We also hnd that p has a tendency towards the upper limit 0.5 at longer integration times. 

ii) We demonstrate t hat the law L oc mP can b e extracted also from compiling data of previous works in the literature 
I Merritt fc Valluri 1 199d ). Kandrup fc Siderii I 2002h ). in galaxies with both smooth and cuspy centers. 

iii) We make a theoretical analysis of the Lyapunov exponents for centrophilic box-like orbits in elliptical galaxies. We 
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Figure 16. Apocentric Poincare surfaces of section with rrihh = 0 at the energies (a) Ej = —204000 and (d) Ej = —195000. (b,e) Same 
as in (a,d), but now with a central black hole of mass = IO^Mq. The xi periodic orbit (c) at the energy Ej = —204000 has a 2:1 
form (two apocentric passages per revolution in the rotating frame), while (f) at the energy Ej = —195000 it has a 4:1 form. 



Figure 17. (a) Tube orbit with energy -216000 around the 2:1 periodic orbit before the insertion of the central mass. The initial 
conditions correspond to a librational invariant curve of the 2:1 island of stability which survives after the insertion of the black hole, 
(b) Same as in (a), but for initial conditions on an invariant curve of the 2:1 island of stability which is destroyed after the insertion of 
the black hole, (c) The orbit with same initial conditions as in (a) remains ordered after the insertion of a central mass = 10® Mq. 
(d) On the contrary, the orbit with same initial conditions as in (b) is converted to chaotic. (e,f) Same as in (a,b) but for two tube orbits 
of the 4:1 resonance. (g,h) The orbits with the same initial conditions as in (e,f) respectively, after the insertion of the central mass 
= 10®Mq. The one of (g) remains regular, while the one of (h) is chaotic. 


demonstrate that the mean Lyapunov exponent can be obtained by independently estimating a) the mean value of the so-called 
‘stretching number’ (=one-step Lyapunov number) at every transit of an orbit from the sphere of influence of the central 
mass, and b) the rate of visits of the orbits to the sphere of influence. In both cases, we find how the various estimates depend 
on m as well as on the orbital energy. Regarding (b), we find two different estimates, according to whether an orbit can be 
characterized as ‘planar’ or ‘3D’. Putting all estimates together, one arrives at a theoretical reproduction of the L oc law. 

iv) In the case of models with central cusps, we find a critical mass scale rUc, such that for central mass parameters 
m < rric the chaotic behavior of the centrophylic orbits is dominated by the central cusp (we call this ‘residual chaos’), while 
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for m > trie an approximate power-law correlation L oc mP is restored, with p « 0.4. The critical mass scale rric as well as the 
‘residual mean Lyapunov exponent’ Lq are increasing functions of the exponent 7 in power-law central cusps p{r) ~ r 

v) We finally explore numerically the correlation between L and m for the centrophylic orbits in disc galaxies with rotating 
bars. In this case, while there can be no box-like centrophilic orbits, we find several quasi-periodic tube orbits around the 
main families of periodic orbits (like 2:1 or 4:1), for which the tube is thick enough so as to pass arbitrarily close to the center. 
These orbits support the rising density profile of the bar at the center. Their initial conditions are close to the last librational 
KAM curve of the islands of stability around their corresponding periodic orbits. Numerically, we observe that only the tube 
orbits around the lowest-order periodic orbits can become centrophylic. Furthermore, for the chaotic counterparts of these 
orbits, after the insertion of the black hole, we numerically recover again a correlation of the form L oc mP. 
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Figure 18. Schematic representation of our model for orbit transits via the sphere of influence of the central mass. An orbit enters 
the sphere at point A and exits at point B. Hyperbolic Keplerian dynamics is assumed in order to estimate the orbit’s local ‘stretching 
number’ arising at the transit (see text). 


Appendix 

We theoretically estimate the local value of the ‘stretching number’ of an orbit transiting the sphere of influence of the 
central mass (Fig llSIl . schematic, as follows: approximating the motion during the transit as Keplerian, in polar co-ordinates 
{r,<j)) with respect to the origin on a plane including the origin and the entry and exit points (A,B respectively), one has 


- = Acos{(j>-(j>o) (39) 

r O'’ 

where C is the local value of the angular momentum (assumed constant during the transit, see section 3). The constant A is 
given by: 


A = 


Gm 


(40) 


where e = -I- 2EC^/mG^, whereas (po is the angle corresponding to the closest approach to the origin, at distance (rmin- 

The orbit has the the same velocity measure vo at the points A and B. We assume that the velocity vector at A forms an 
angle u with the horizontal axis. 

With the above conventions, to a given entry angle ^ corresponds a given exit angle (j>' from the sphere of influence. After 
some algebra (taking into account Eg. (1391) as well as the preservation of energy and angular momentum), we find: 


= 2 cos ^[ 


Vo sin (0 -I- u) 


V(G Tfilvm'EQ sin (0 -|- ti)) -|- Vq — ‘ZGTnjVri 

-■ 

TmVQ sm (0 -I- U) 


(41) 


The local value of the stretching number can now be estimated by the difference in 
slightly different angles (f>. Taking the derivative of eq. m we have: 


for two nearby orbits entering at 




Vo sin (0 -b u) 


V(G n?'//' 7 Ti,no sin ^G^TnjT'n 

( 1 - 


rmVQ sin^ {cj>-\-u)' 


(42) 


Re-orienting the frame of reference, without loss of generality the parameter u can be set equal to zero. The quantity /d(j)\ 
is the measure of the stretching number a for the transit, motion inside the sphere One can check that for typical energies 
of centrophilic orbits, one has 2mr-mVo << rfvo- Then, Ea. ()42l) reduces to: 


d(j )' . M _ 1 -b mcsc^ (j>/rmVo 

d(j) 1 -b m? csc^ (j)/r'^VQ 


( 43 ) 
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The mean stretching number of a transit at given velocity vo, with respect to all possible angles (p can be estimated by the 
integral of the quantity log \ d(f)'/d(j)\ over all possible angles 


log 




dcj) 
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TT 
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TT 
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TT 


'r/2 ,il 


dcj) 


^7r/2 


log|l-2 


1 H-^ csc^ cj) 


1 + ■ 


-\d(f) = 


r-^/2 


, I/-I o“0 

log 1(1 - 2 3 2 , )\d<l> = 

Og + sm (f> 


(44) 


^7r/2 


log (2 


Qq + sin^ (j> 


, 1 + ao 


2 , „;„2 ^ - l)rf<(> = log (2 t^ - 1) + 2 log 




Qg + sm 


1 + Qg 


Q?0 


+ \/l + “§ 


where qq = m/r^Ug. At mass ranges 10“® < m < 10“^ the value of ao is in general a small quantity, except for small energies 
(E < 0.03 in our unit s), which, however, can b e readily checked to correspond to orbits residing always within the sphere 
r = r„i, i.e. pyramids ( Merritt fc VasilievI ( 20111) '). Excluding such orbits, we set ag = 0 in the previous equation and obtain: 


log 


dcj)' 


dcj) 


~ log(l + 2ao) + 2 log 


(1 + V^^g^) 


Q) \ 


= log(l + 2ao) + 2 log 


(1 + 


0/ “h 1 


(45) 


Expanding Eg. (145 II with respect we find: 


log I —^ I — 2ao + 2 log (1 + yaoV^) — 2ag -j- ao 

dcj) 


-m + 


V22 


(rmV^) (r^ug)i/2 


^/m ~ 


(rmWg 


-%/m . 


(46) 


In Ea. (l46ll the quantity depends on m (r^ oc Furthermore in the limit of the sphere of black hole’s influence, energy 

is approximately equal with the kinetic energy (vq oc E), whereby the last expression can be written in the form: 

_ dcj)' 

i.e. we arrive at the estimate of Eq. nil). 


( 47 ) 





































